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Abstract 

In this paper we consider the effect of self-localization of a quantum impurity in Bose-Einstein 
condensate. Space correlation function of the impurity is evaluated with the help of the imaginary- 
time path integral approach. Employing the Feynman's variational method we calculate the im- 
purity correlation function as well as the energy of the system associated with the impurity. The 
effect of self-localization predicted before within Gross-Pitaevskii approach is recovered by our 
treatment. The strong coupling regime with negative ground state energy is reached by variational 
method, and corresponding correlation function is calculated. 
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I. INTRODUCTION 



A theory of atomic size impurity in Bose Einstein fluids has become a subject of intensive 
research in the last decade. The theoretical development is largely motivated by recent 
experiments. MIT group pQ has studied the dynamics of impurity atoms in Bose Einstein 
condensate (BEC) of ultracold atoms in magnetic traps. Another set of experiments studies 
the microscopic superfluidity of liquid helium by spectroscopic mesurements on molecules 
imbedded in helium droplets [21 13]. It is generally observed that microscopic impurity particle 
interacting with Bose liquid/gas behaves as a free particle with effective mass increased by 
its interaction with Bose fluid. Such dissipationless quantum motion is observed both for 
translational motion of particles in BEC and for rotations of molecules in superfluid helium. 

A question that received attention of theorists has been the structure of the ground state 
of impurity in BEC, in particular, the possibility of the so called self- localization of an 
impurity: the appearance of a bound state of an impurity with BEC despite the purely 
repulsive interaction potential. Such behavior has been investigated by a number of authors 
using linearized Gross-Pitaevskii equations |U El El E|- This method allows to obtain the 
non-linear imaginary time Schrodinger equation for the impurity particle. The energy of 
the particle becomes negative above certain value of coupling strength that depends on the 
BEC/particle and BEC/BEC interaction potentials. Extension of these results within fully 
quantum treatment is an open question and is the main subject of present work. 

In this paper the matrix element of the reduced density operator of impurity is calcu- 
lated which allows to compute the energy and the space correlation function of impurity. 
The system of our interest is the particle interacting with the gas of uncoupled Bogoliubov's 
excitations. We formulate this matrix element via imaginary time path integrals. The inte- 
gral over BEC trajectories is explicitly eliminated and the non-Gaussian functional integral 
over impurity trajectories is left to be calculated. We treat this integral using the varia- 
tional approach developed by Feynman for the polarons in polar crystals [8]. This method 
was successfully applied to the problems described by the Frohlich-type Hamiltonians such 
as polarons, electron-plasmon interaction [9], nucleon-meson interaction [TO], many body 
fermion problems [TT] . and many others. Since the Hamiltonian of our model (one par- 
ticle interacting with uncoupled Bosonic modes) is of the Frohlich type, it lets us believe 
that such treatment is successful in our case. We obtain the free energy of an impurity as 
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a function of coupling constant. The energy is increasing in the weak-coupling regime in 
agreement with the regular perturbation theory treatments. With further increase of cou- 
pling it reaches maximum and decreases becoming negative in the case of strong-coupling. 
This indicates the existence of the bound state. By computing the correlation function we 
obtain the localization radius as a function of the coupling constant. We show that the 
critical localization radius is of the order of magnitude of the inverse critical momentum 
above which the dissipation of particle motion takes place in real time dynamics. 

The Plank constant % and the Boltzmann constant ks are set to unity throughout the 
paper. 



II. MODEL HAMILTON! AN AND STATEMENT OF A PROBLEM 

As it has been mentioned in the Sec. [T| we are going to concentrate on the case of dilute 
Bose gas at the temperature mach less than a temperature of condensation. So the gas 
of weakly interacting Bose particles with mass m can be described as a gas of uncoupled 
Bogoluibov's excitations [121 H3], i-e. the Hamiltonian of BEC has diagonal form and reads 

*- = I»*4, «w = ^(^+w). (i) 

Here bosonic operators and create and annihilate the collective excitation with 
momentum k and with the Bogoliubov's spectrum e(k) which has the phonon-like behavior 
at low momenta, i.e. eik — > 0) = kc, where c is the speed of sound. 

Next, we consider a quantum particle with mass M interacting with the BEC discussed 
above. The Hamiltonian of the whole system has the form [H] 

H = gn + H B + J2 + Yl 7k (K-^B£ + a+ +k a q fl k ) , (2) 

q q,k^0 

g I Nk 2 
7k ~ V\J 2me(k) ' 

The bosonic operators a q and a q in the above expression create and annihilate a single 
particle with the mass M in the state with momentum p. The last term in Eq. ^ describes 
the particle-BEC interaction with the coupling constant 7 q depending on whole number of 
particles in bose gas N in the volume V = N/uq. Also, the coupling strength depends on 
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the coupling constant, g, denned as a zero Fourier component of inpurity/BEC interaction 



potential. This constant is related to the speed of sound as c = \/gn /m B . The first term 
in ([5]) represents the first order correction to the energy due to the particle-BEC interaction. 
In the interaction part of the Hamiltonian as well as in the Hamiltonian of the free BEC 
we neglected the terms responsible for the interaction between the Bogoliubov's excitations. 
This approximation remains valid if the single impurity alters the surrounding BEC only 
slightly which is always the case in the macroscopic limit. 

In order to describe the statistical properties of the impurity that stays in thermal equi- 
librium with the BEC environment, we will calculate the space correlation function defined 
as a matrix element of reduced density operator 

1 



C(x,x') = -(x|Tr B e-^|x'> • 



Z 



(3) 



Here j3 is the inverse temperature, the Hamiltonian H is defined by Eq. (|2j), the trace is 
performed over all states of BEC and Z is the partition function. Due to the momentum 
conservation in the whole system the reduced density operator is diagonal in momentum 
space. Thus 



C(x 



ZV 



(4) 



p p 

So our task is to calculate the correlation function as the reduced density matrix element 

between two states of the impurity with the momentum p 



1 



C P = p p ((3) = -Tr B (p|e-^ip) = -Tr B (0|a p e-^a+|0) 



(5) 



The above matrix element can be represented as the imaginary-time coherent-state func- 
tional integral [15] 



cm 



1 

z 



£>[K(t)},k(t)}] / n 



dbldbk 
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D[{bi(r)},{b k (r)}}a p (pym 



x oxp - Y,\b*\ 2 + Y, b UP)bk((3)-S 

k k 

with the imaginary-time action of the whole system 



(6) 



S 



dr 



£ M r KM + ^(q)<(r)a q (r)) + £ (b k (r)b*( 

k 

« q (r)6 k (r)) 



<k)^(r)fe k (r)) 



+ E 7k ( a q+k( r ) a q( r ) 6k ( r 



«q-kl r J 



(7) 



The integral ^ must be evaluated with the following boundary conditions for the trajecto- 
ries 

6k(0) = & k , b*M = bi , (8) 
a q (0) = a*(/3) = . (9) 

Because of the linear dependence of the interaction part of the action on the b^ir), b^(r) 
trajectories the BEC degrees of freedom can be integrated out immediately using the formula 

^D[{b* k (r)},{b k (r)}]e W - £ jf dr{e(k)b* k (T)b k (r) - j k (r)b k (r) -j k (r)6*(r 

= Z B exp [ dr J* dr'T^r - T')j k (r)j* k (r')\ . (10) 

The function i\(t) represents the imaginary-time propagator of free Bose field with the 
spectrum e(k) 

T k (r) = (0(r) + n k )e-"« , n k = ^ _ ^ , (11) 

where function O(r) denotes the Heaviside step function, and Z B = Yl^n k e^ £ ^ is its par- 
tition function. 

Thus, after tracing out the BEC we obtain the functional integral representation of the 
reduced density matrix of the relevant particle that includes the integration over the particle 
trajectories only 

PM = ^ J D[i<(r)}, {a q (r)}]a p (/3)a;(0) exp [S P + Sj] (12) 
Here we defined the free imaginary-time action of relevant particle 

S P = ! drJ2 MtK(t) + £(q)a q (r)a q (r)) (13) 

JO 



The integral (12) is of the non-Gaussian type and its non-Gaussian part Sj is due to the 



impurity-BEC coupling and reads 



rp rp 

Si= dr dr' V 7k a q+kW a q( r ) r k(r-r> q ,_ k (r> q /(r') . (14) 

Having the action of the impurity, we can investigate its classical dynamics described 
by corresponding Euler's equations of motion. For this purpose let us write the action 
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in the coordinate representation by changing the variables of the functional integration as 
ak(r) = ^= J dre~ t]ir ilj(r,T). Now the action of the impurity surrounded by BEC has the 
form 



S imp = S P -Sj = J J drr(r, r) (J- - ^ (r , r ) _ (15) 

J dr J dr' J dr J dr'ip*{r, r)^(r, r)A(r — r', r — r')tp*(r', t')i/)(t', t') , 



A(r - r\ r - t') = Ap J dk^r - r')e jkr (16) 

Varying the above action one gets the following equations of motion 

^(r,r) = ^^(r,r)+ / dr' J dr'ip(r, r)A(r - r, r - t')^*^', r')^(r', r) , (17) 

-V>*(r,r) = A^( r ,r) + J dr' J drty*(r, r)A(r - r\ r - r')^(r', r>(r', r') , 

with symmetrized kernel 

A(r — r',r — t) = A(r — r', t — r') + A(r' — r, t — r) . 

We will seek the stationary solutions in the form if)(r,r) = ip s (r)e~ TE and ^*(r,r) = 
ip*(r)e TE . So the stationary equation for the impurity wave function is 

E^ s {r) + ^s(r) + J C /rV s (r)^(r')^(r / )A s (r - r') = , (18) 



where 



A s (r -r') = J dr'A(r -r',r- r') . (19) 



Using the definition of the BEC-propagator Eq. (11) together with the Eq. (16) one 
obtains the non-linear kernel A s , given by 

= (20) 

7r r 



The equations (18) with the kernel (20) represent the particle self- interacting through the 
screened Coulomb potential. The same equations for the impurity in weak-coupled BEC can 
also be obtained directly with the help of linearized Gross-Pitaevskii equations with more 
restricting initial assumptions [5] such as Hartree approximation for the impurity/BEC wave 



function and using the real-valued BEC wave function. In Ref. [5] authors investigated the 



lowest energy solution of the above equations. They found that the equations (18) have 
the localized solution with negative energy in the strong coupling regime, i.e. when the 
impurity/BEC interaction is strong enough to compensate the high kinetic energy of the 
localized state of the impurity. 

The aim of this paper is to construct fully quantum description of the model system based 
on the calculation of the imaginary-time quantum propagator. The quantum formulation 
of the problem without the restriction of a mean-field approximation will be given in the 
following section. 



III. PATH INTEGRAL FORMULATION OF THE CORRELATION FUNCTION 



Our goal is to describe the effect of self-localization from purely quantum point of view. 
For this purpose let us return to the functional integral for the correlation function with 
the full impurity-plus-BEC action of the form Our aim is to represent the correlation 
function as a single Feynman's path integral over the impurity trajectories. As in the 
previous section we will use the coordinate representation for impurity trajectories. So the 
correlation function is given by 

C(x-x',/3) = ^D[{bt(r)},{b k (r)}] J D[{^(r, r)}, {V>(r, r)}]^(x', /3)V*(x, 0)e~ s (21) 

where the action is 

S = ! drJ2(k(r)b*(T) + e(k)bUr)b k (r))+ I dr f dr L\r, r)^(r, r) - ^*(r, r) ^(r, r 

x jf dr J ^(r,r)^(r,r)^7 k (6 k (r)e ikr + 6^(r)e- jkr ) . (22) 
Let us introduce auxiliary external sources into the action and define the following functional 



F[f(r,r),j(r,r)} = J D[{^*(r, r)}, {^(r, r)}] 



(23) 



x exp 



- / dr I dr [^*(t,t)K(t,t)iI>(t,t) - j(r,r)^(r,r) - j*(r, r)^*(r, r) 



K(r,r) 







dr 2M 



A £ 7k (6k(r)e ikr + 6 k (r) e - kr ) = A + H(r, r) . 

k 



(24) 
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The functional integral in the right side of Eq. ( 23 ) has to be evaluated with the boundary 



conditions ip (r,0) = ip*(r,f3) = 0. Since the integral in (23) is Gaussian, it can be done by 



the stationary phase method. Thus for the stationary trajectories one finds 

V>(r, r) - H{r, r)^(r, r) - f (r, r) = , 
-ip*(r,T) -H(r,T)r(r,r) - j(r,r) = . 



(25) 



The formal solution of the above equations can be written in the following form 

^(r,r)= / dr'exp( [ H(r, s)ds ) j*(r, r') . (26) 



Finally, substituting this solution into the integrand in (|23|), for the functional F we have 

(27) 



F[f(r,T),j(r,r)} = exp ^ dr ^ drj(r, r)^ fl (r, r)^ 



rP rP 

vxy>[ l dr J dr'Q{T -T')(j{v,T)\e~^ H{ ^ s)ds \f{r,r')) 

Now we can note that the correlation function can be written with the help of the generating 
functional, namely 

5 2 



C(x-x',/3) 



F[f(r,r),j(r,r)] 



(28) 



j=j*=0 

where the averaging is performed as the integration over BEC trajectories with the free 
action of BEC as 



(■■■) = ^4>D[{h k ir)\.{b k (T)}] ... oxp 



P 



dr^(6 k (T)fo*(T) + e(k)&*(r)6 k (T)) 



• (29) 



Calculating the functional derivative in Eq. (28) and writing one particle propagator 



(j(r,r)|e /t' w ( r ' s ) ds |j*(r, r')) as Feynman's path integral, for the correlation function we 
have 

/ /-x(/3)=x' \ 

C(x-x',/3) = ( / D[x(r)]e-M , (30) 

\7x(0)=x / 

where the single particle action depends on the coordinate trajectory x(r) 



S ,, 



dr 



Mx 2 



+ $>k (6 k (r)e ikx ^) + b* k (r)e~^) 



(31) 
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Now using the formula (10) we can eliminate the integration over BEC trajectories from the 



equation (30) and get the correlation function as a single Feynman's path integral for the 
impurity 

,x(/3)=x' 



C(x-x',/3) 



D[x(r)]e- 5fl 



(32) 



x(0)=x 



where the impurity action reads 



S f 



/ dr—- dr rfr / V 7 2 r k (r-r')e 
Jo 1 Jo Jo ,, 



,/3 r /3 



ik(x(r)-x(r')) 



(33) 



This form of the single particle functional integral is considered in the following section. 



IV. FEYNMAN'S VARIATIONAL APPROACH 



The imaginary-time correlation function in the form of path integral (32) provides the 



exact description of the statistics of the impurity surrounded by the degenerate BEC at low 
temperature. In this section we will calculate this integral approximately using Feynman's 
variational approach to the polaron problem [S] . The idea of the original method is to replace 



the real action of the impurity (|33|) with the trial action of the form 

- p , Mx 2 , 

dr— h dr 

'0 2 J Q J Q 



r p Mx 2 r p r p 

S T = / dr— + / dr / dr'Q(r-r')(^(r) 
Jo 2 J Q J Q 



xir 



')) : 



(34) 



where 



Q( T - T >) = g (6(r - r') + nJ)e~^ T -^ 



(35) 



It is also useful to write the trial action in the extended form as 



e~ Sr ~ 



j D[y{r)}e- SET 



S 



ET 



dr 



Mx 2 (r) my 2 (r) mu 2 (y(r) - x(t))' 
2 2 + 2 



(36) 
(37) 



Performing the integration over trajectories y(r) one recovers the trial action in the form of 
Eq. (34) with Q — rnoo 3 /4. Thus in this variational treatment we replace the interaction 



with the original BEC environment by the interaction with single trial particle with the 
mass m. 
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The trial action of the form (34) is assumed to be the zeroth order contribution while 



the difference between real action (33) and trial one has to be considered as a perturbation. 



Thus one has to write the following expansion of the correlation function 

(1 + S T - S r ) 



C t (x-x',/3) 
where we have defined two kinds of averages 

rx(/3)=x' 



(1 + St — Sr)q 



(38) 



(...) 



D[x(r)]... e 



-St 



(-)( 



x(0)=x 



x(/3)=0 



£>[x(t)]... e~ ST 



(39) 



x(0)=0 



The trial correlation function Ct in Eq. (38) still depends on two variational parameters 



Qo and u (or m and to). Thus, in accordance with the principle of minimal sensitivity 
one has to minimize this expansion with respect to the variational parameters, and the 



extremum point will give the best variational approximation for the path integral (32). In 



order to proceed with the expansion (38) we will have to calculate four functions, namely 
Si 



and 



S? 



yO 



[ dr [ dr'g(r-r')<(x(r)-x(r')) 2 > , 
Jo Jo 

I dr ! dr' V T k (r - r')< exp [zk(x(r) - x(r'))] > 
Jo Jo k 

fdr rrfr'g(r-r')<(x(r)-x(r')) 2 > , 
Jo Jo 

fdr fdr'^Mr ~ r')(exp [ik(x(r) - x(r'))] > 
Jo Jo ,, 



(40) 



(41) 



The expansion (38) in terms of functions E becomes 

(1) + S 1 + S 2 



Cn 



(l)o + S? + SQ 

>2\ 



(42) 



In order to calculate the average ((x(r) — x(r')) ) and (exp [ik(x(r) — x(r'))] ) that 
enters the functions Si and S 2 , respectively, one can use the generating functional of the 
form 



/•x(/3)=x' / p/3 

W[j(r)] = j> D[y(r)} / D[x(r)] exp -S ET + / dsx(s)j(s] 

^x(0)=x V Jo 



(43) 



So the second average (exp [zk(x(r) — x(r'))]) can be evaluated by setting j(s) = zk(5(s — 
r) — 5(s — r')) while the first one is obtained by differentiating the second one twice over k 
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at k = 0. The functional integral in (|43j) can be evaluated by introducing the new variables 
q(r) = x(r) — y(r) and r(r) = (Mx(r) + my(T)) /{M + m). Substituting this replacement 



into (36) one gets the action Set in the form 



S 



ET 



dr 



mM 



(M + m)r (r) //q (r) mw q (r 



+ 



+ 



j(r)r(r)--^j(r)q(r) 



(44) 



m + M 

The above action now describes harmonic oscillator and free particle in presence of the 
external source and the corresponding functional integral can be easily evaluated. Here we 



will write down the result of the integration in (|43|) 
H/[j(r)](x,x') = W exp (-(x - x') 2 



Mu Qu , fi/3 
+ — coth — 
2/3m 4 2 



(45) 



+ / drj(r) 



o 



' it t u a a sinh fir- sinh fi(/3-r) , 
(X " X W + ™ X+ ^ (X + X) + 2M simTfi^ (X " X) 



+ / <*r / dr'j(r)j(r') 
Jo Jo 



r ' /9_r at >\j_ T & ~ T ' at ' \ 

-U{T — T ) + V{T — T) 



2(3 M + m 



2(3 M + m 



11 If sinh - r) + sinh fir + sinh fi(/3 - r') + sinh fir' ^ 



4fiM 2 tanh ^ V 



sinh 0,(3 



H coshO(/3/2- |r-r'|) 
_ 4fiM 2 sinh^ 



(46) 



where we have introduced the frequency fi = {jjyfmfji. The prefactor W is the constant 
coming from the integration over the deviations around the stationary trajectories, so it is 
independent on the end-points x and x'. 

So for the second average (exp [zk(x(r) — x(r')jl) one gets 

"Mu 2 Qfx , np~ 



(exp [ik(x(r) — x(r'))]) = lVoexp< — (x — x') 



2/5^ + X C ° th 2 



(47) 



+zk(x' — x) 



to 2 t - r' 1 fi 2 - uj 2 sinh fir - sinh fi(/3 - r) - sinh fir' + sinh Q((3 - r') 

fi 2 p ' 



2 fi 2 



-k 2 



1 U) 2 I \T- 7 

r — r 1 — 



2Mfi 2 ' 







sinh fi/? 



1 fi 2 - to 2 cosh ^ - cosh fi (p/2 



r — r 



'I) 



2Mfi fi 2 



sinh 
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Below we will be always interested in low temperature limit (3 — > oo. Calculating the first 
average ((x(r) — x(r')) 2 ) with the help of the above expression, for the function S x one gets 



Si 



: l\2 , 

x-x + 



cxp 



X — X 

4 v ; 



(48) 



2u;0 3 v " ' MwO_ 

The calculation of S 2 can not be done in closed form and requires some numerical calculations 
that will be performed in the next sections. 



V. ENERGY OF THE IMPURITY 



In this section we will investigate the ground state energy of the impurity as a function 
of the impurity-BEC coupling strength. Let us define the free energy of the system E as 

rx(/3)=0 



-f3E 



D[k(t)} e" SR = (e ST ~ SR ) . 



(49) 



x(0)=0 



In the zero temperature limit free energy coincides with the ordinary energy and shows the 
difference between coupled and uncoupled impurity/BEC system. 



The expression (49) up to first cumulant reads 

(St — Sr) 



-PE(n,u) 



(l)oexp 



(l)c 



E(Q, u) = E + E 1 + E 2 



(50) 



The zero energy contribution Eq is 

,x(/3)=0 



-PE 



D[x(r)]e- 



W[j = 0,x-x' = 0] = iy 



(51) 



'x(0)=0 

and the first order contribution terms are defined as 



1 



S? , E 2 



PW 



-eS 



2 • 



(52) 



where S" 2 is defined by Eq. (41). While the terms E\ and E 2 can be evaluated directly 



using the generating functional, the first term Eq can not be calculated explicitly from the 



formula (51). In order to find E one has to note that 

fdr f dr'Q(r-r'){(^(r)-^)) 2 ) i 
Jo Jo 



d 

-m— — e 
dm 



-f3E 



Using the result (48) we can get following equation for zero order energy 

3 uo 2 



dE n 



dm 4Mfi 



E (m = 0) = 0. 



(53) 



(54) 
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The solution of the above equation is 

3 



(55) 



Now it is left to calculate the contribution E 2 . With the help of the formula (47), in the 
zero temperature limit one can give the following expression for E 2 



Eo = — 



k 2 fu 2 T 



- e(k)r -yl^M + 



n 2 



(56) 



First, let us consider the weak coupling limit, i.e. Qo — > or m — > so that Q ~ uj. Thus 
we have in the lowest order Eq + E\ = and the energy E is given by the first expansion 
term of the ordinary perturbation theory 

1 



E 



weak 



E. 



weak 



k C ( k ) + 2M 



(57) 



In order to prevent the ultraviolet divergence in the above sum over momenta one has to 
renormalize the coupling constant g according to the second order Born approximation. 
For this purpose we must add the zero order energy term gn from the Hamiltonian ^ to 
the energy expression and expand the constant g in powers of scattering length a for the 
impurity-BEC interaction 



2%a , 
9 = II 



2a 



dk 



m B M 



(58) 



id,. \ it J I rriB + M 

The weak coupling energy reexpanded in powers of scattering length in now finite and reads 



E weak = m B c 2 a(l + z) 



dl 



1 



l + z 



(59) 



7 1 + 4// 2 (yi+4/z 2 + ^ 

Here we have introduced dimensionless coupling parameter a = 4an/(mBc) and the mass 
ratio z = m B jM. Below we will always use the renormalized expression for E 2 which is 

4a 2 n 



E 2 + 



m r 



dk 



(60) 



Finally, the whole expression for the ground state energy can be written in the integral form 
3{Q-uj) 2 



E ren (Q,uj) 



4 n 

+ a(l + z) 



(61) 



dl 



1 - 



l + z 



2 y/l+4/l 2 



ds exp 



l 2 <Tl(s) 
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where 

Ms) = s^^+sz^ + ^^^(l-e-^) , (62) 



The integration in (61) is performed over dimensionless variables s and /, and the energy 
E(fl,oo) as well as the frequencies Q and lo are now measured in the units mgc 2 . 

Following the principal of minimal sensitivity [TB] (or Feynman- Jensen variational in- 
equality in case of the energy) we seek the extremum point of the energy with respect to the 
variational parameters VL and lo. Fig. shows the energy as a function of f2, lo for the mass 
ratio z = 4/3, which corresponds He3 impurity in He4 BEC, and for the value of coupling 
constant a = 5. This dependence has the form of gutter having a minimum with respect 



to Vl and not having it with respect to lo. As it can be seen from the Eq. (61), the strong 
coupling asymptotic of the energy does not depend on lo and behaves like E strong ~ —a 2 
requiring the optimization with respect to Q only. On the other hand, the weak coupling 
regime is independent of both variational parameters. Thus only the intermediate coupling 
region is sensitive to the choice of the second parameter lo. On Fig. ^ we plotted the energy 
as a function of the coupling constant a optimized with respect to Q for different values of 
the second variational parameter lo. It turns out that at some value of the coupling param- 
eter a the energy becomes negative indicating the existence of the bound impurity/BEC 
state in spite of the repulsive interaction potential. One can see that the optimized value of 
the energy rapidly converges to some minimum with increasing lo becoming insensitive to lo 
at the values of lo > 50. Thus the use of this value is safe within variational treatment. So 
the critical value of the coupling constant a c defined as E(a c ) = decreases from a c ~ 3.2 
at lo ~ 1 down to reliable value a c ~ 1.8 at high lo. Next, Fig. ^ shows the critical value 
of the coupling constant function of the mass ratio m^/M. As one would expect 

the impurity never self-localizes if its mass is much bigger than the one of the Bose particle 
and binds with BEC in the weak coupling regime in the opposite case of the small impurity 
mass. 
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VI. IMPURITY CORRELATION FUNCTION 



Now let us directly consider the correlation function C{x). Using expression for the 



function Si, Eq. (|48|), one can rewrite the expansion (38) up to first cumulant as 

Mfl ftt 2 -u 2 ^ 2 



CV(x,/?^oo) 



cxp 



(3E 2 



x 2 + S2(X) 



(63) 



tt 2 J " ' W (x) 

The above function is normalized in such a way that Cr(x = 0) = 1. Then, with the help of 



the Eqs. (45, 47) the correlation function can be written down in the integral form, namely 

exp 



CV(x, (3 — > oo) 



x exp 



n 2 



PaAs-s') 



2 [I + Z) 

x + a — 



P 



ds / ds' 



r 2 



o v/l + 4/Z 2 



sin(ix(72(s — s')) 

lx<Jo(s — S') 



where the function <7i(s — s') is defined by the Eq. (62) and 



J 2 s 



lfi 2 

2 fi 2 sinhfi/3 



ft 2 /?' 

x (sinh fis — sinh Q(/3 — s) — sinh Qs' + sinh — s 1 )) . 



(64) 



(65) 



As in Eq. (61), the integration on the right hand side of the above expression is performed 



over dimensionless variables s, s', and I, temperature 1/13 and freqiencies Q and uj are 
measured in units m#c 2 while the variable x is defined in units 1/mgc. Finally, in order 



to find the best approximation for the correlation function the expansion (64) has to be 



optimized with respect to the pair of the variational parameters Q, uj. In general, if the 
optimization procedure is directly applied to the expansion for the correlation function, the 
optimized values of variational parameters should depend on the coordinate x. However, it 
turns out that instead of the variational optimization of the coordinate-dependent correlation 
function one can use the optimal values of the variational parameters obtained from the 
optimization of energy. These coordinate-independent values then has to be substituted 
into the expansion for the correlation function giving the variationally improved result, 
i.e. convergent strong coupling expansion for the correlation function. This scheme of 
the variational perturbation theory was successfully applied for the variational solutions of 
time-dependent equations in non-linear dynamics [T7] [IB"]. 

We will also define the mean value of x 2 representing the square of the radius of localiza- 
tion or polaron radius as 

2n J x 2 C(x)<ix 



(x 2 ) 



/C(x)dx 



(66) 
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Fig. (|4]) shows the mean square of the radius of localization as a function of coupling 
constant. It has inverse dependence on coupling strength. The inset shows the correla- 
tion function itself for different values of coupling strength. The correlation function has 
Gaussian-like shape regardless of a which means that formally the particle is always local- 
ized. This is the consequence of the choice of the trial action used in this method. How- 
ever, this correlation function should converge to the true-one at large values of interaction 
strength. It is interesting to examine the localization radius at the coupling constant a e at 
which the energy has maximum, i.e. dE(a)/da\ ae = 0. We notice that (x 2 (a e )) ~ 1/m^c . 
This can be expressed through a so called healing length of the BEC £ = l/zp c , where p c 
is the critical momentum of the impurity above which the dissipation takes place in case of 
real-time dynamics in accordance to Landau's criterion. This leads to a/ (x 2 (a e )} ~ l/p c , or 
using Heisenberg uncertainty for the localized particle \J (p 2 ) ~ p Cl i.e. the average momen- 
tum of a particle at the point of self localization is of the order of the critical momentum. 
One can make physical sense of this fact by noticing that this means that the tendency of 
self-localization appears when the exchange of energy between particle and BEC becomes 
possible. 



VII. CONCLUSION 



In this work we consider a problem of self- localization of impurity in degenerate BEC. 
This problem has been considered previously using the mean- field approach. Here we use a 
full quantum description of the ground state of impurity surrounded by BEC. The variational 
perturbation method is employed to calculate the imaginary-time propagator of impurity in 
degenerate BEC. The free energy and the spatial correlation of impurity in BEC is obtained 
as a function of coupling strength. Our results point to a possible relation between self 
localization and the real-time particle/BEC energy exchange. In this work we explicitly 
used the degeneracy of the BEC. However, the methodology developed in this paper can be 
extended to Bose systems without an assumption of diluteness. The latter can be done by 
using the four-point Green's function of the strong coupled Bose liquid [12] in the leading 
expansion term instead of degenerate BEC propagator used in this work. 
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FIG. 2: Energy optimized with respect to Q as a function of the coupling constant a for different 
values of second variational parameter u. 
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FIG. 3: Critical value of coupling constant a c as a function of the mass ratio ms/M. 
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FIG. 4: Mean value of x 2 as a function of coupling constant, x is shown in units l/m^c. Inset - 
correlation function c(x) for different values of coupling constant a = 0.5 (upper curve), a = 1,5 
and a = 10 (lower curve). 
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